Volume terms for charged colloids: a grand-canonical treatment 
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We present a study of thermodynamic properties of suspensions of charged colloids on the basis 
of linear Poisson-Boltzmann theory. We calculate the effective Hamiltonian of the colloids by 
integrating out the ionic degrees of freedom grand-canonically. This procedure not only yields 
the well-known pairwise screened- Coulomb interaction between the colloids, but also additional 
volume terms which affect the phase behavior and the thermodynamic properties such as the osmotic 
pressure. These calculations are greatly facilitated by the grand-canonical character of our treatment 
of the ions, and allow for relatively fast computations compared to earlier studies in the canonical 
ensemble. Moreover, the present derivation of the volume terms are relatively simple, make a direct 
connection with Donnan equilibrium, yield an explicit expression for the effective screening constant, 
and allow for extensions to include, for instance, nonlinear effects. 



I. INTRODUCTION 

Colloidal suspensions are multi-component systems that 
consist of mesoscopic colloidal particles dispersed in a 
molecular solvent. Often other chemical components are 
present as well, e.g. ions, polymers, or proteins. Pre- 
dicting or understanding the properties of such mixtures 
from a microscopic perspective is generally complicated, 
because the large asymmetry in size and charge between 
the colloids and the other components in practice inhibits 
a treatment of all the components on an equal footing. 
The standard way out is to regard the suspension as an 
effective colloids-only system, in which all microscopic 
degrees of freedom of the "medium" (solvent, ions, poly- 
mers, etc.) are suitably averaged out. 

For instance, in the case of colloidal hard spheres 
in a medium with non-adsorbing ideal polymers (ra- 
dius of gyration R g ) the so-called depletion effect Q 
leads to effective attractions between pairs of colloids 
at surface-surface separations less than 2R S , and in the 
case of charged colloidal spheres in an electrolyte with 
Debye length re -1 the effective interaction between a 
colloidal pair at center-to-center separation r is gener- 
ally written as a repulsive screened-Coulomb potential 
oc exp(— nr)/r. The advantage of such a one-component 
viewpoint is that all of the machinery of classical one- 
component fluids (integral equations, perturbation the- 
ory, simulation, etc.) can be employed to study the prop- 
erties of colloidal suspensions, but only after a reliable av- 
eraging over the medium has been performed. Perform- 
ing this averaging explicitly is generally a tremendous 
statistical mechanics problem, that can only be solved 
approximately in most cases 0- 13-01 • 

One important problem is that the effective colloidal 
interactions are not necessarily pairwise additive, i.e. 
triplet or higher-order many-body potentials may ap- 
pear even if the underlying interactions in and with the 



medium are strictly pairwise. On physical grounds one 
generally expects the breakdown of pairwise additivity 
of the effective interactions if the typical length scale in 
the medium is of the order of the typical colloidal length 
scale, e.g. the colloidal radius a. For colloid-polymer 
mixtures it was indeed shown that equal-sized colloids 
and polymers (i? g = a) have bulk and interfacial prop- 
erties that differ dramatically from pairwise predictions 
and charged colloids in an electrolyte were shown 
to exhibit non-negligible effective triplet attractions on 
top of the pairwise repulsions 8\ at (extremely) low salt 
concentrations where k" 1 ~ a. 

In this paper we will focus on a description of effec- 
tive interactions (or the effective Hamiltonian) in bulk 
suspensions of charged colloids. The classical theory for 
these systems dates back to the 1940's, when Derjaguin 
and Landau 9] and Verweij and Overbeek ^3] indepen- 
dently calculated that the effective potential between two 
identical homogeneously charged colloidal spheres (ra- 
dius a, total charge — Ze with e the proton charge) in a 
bulk medium with dielectric constant e and Debye length 



is given by 



3C 



V 2 (r) 



Z 2 e 2 /exp(fta)\ exp(— nr) 



1 



r <2a 



r > 2a. 



(1) 



Here and in the remainder of this paper we ignore the 
dispersion forces, and we recall that K is defined as 



K = 8it\bc s 



(2) 



in the case of a 1:1 electrolyte with total ion concentration 
2c s far from the colloids, where Ab = /5e 2 /e is the Bjer- 
rum length, (3 = and T the temperature |allf|. 

It is well established by now that many properties of 
suspensions of N 3> 1 charged colloids in a solvent vol- 
ume V (density n — N/V) can be understood on the 
basis of the pairwise effective Hamiltonian 
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(3) 
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where denotes the position of colloid i — 1, . . . , N, 
and where i?y = |Rj — Rj|. For instance, the thermo- 
dynamic equilibrium properties and the phase behavior 
follow from the Helmholtz free energy F% (N, V, T) , de- 
fined as the classical canonical phase space integral 



exp(-/3F 2 ) 



1 



dR N exp(-f3H 2 



Etr c exp(-/3£f 2 ), 



(4) 



where V is an irrelevant constant volume (accounting for 
the internal partition function of the colloids) that we in- 
clude for dimensional reasons, and where we introduced 
the short-hand notation tr c for the classical canonical 
trace over the colloid degrees of freedom. On the ba- 
sis of Eqs. (QJ, 0, and Q) one can explain many ex- 
perimental observations, including the crystallization of 
(essentially) hard spheres (Z = or kq > 1) at packing 
fractions rj = 4ira 3 n/3 > 0.5 into an FCC crystal (TlHy 
the crystallization into BCC crystals for sufficiently soft 
spheres (l^. H^|. the measured osmotic equation of state 
pjj 03, structure factor llfllial, ra dial distribution func- 
tion pair interactions [2CLl2lj. and many other col- 
loidal phenomena. It is therefore fair to state that the 
DLVO theory as described by the Eqs. (O, and <@J 
is one of the corner stones of colloid science. 

It is, however, also fair to add that not all experi- 
mental observations are in (qualitative) agreement with 
DLVO theory. For instance, the experimental observa- 
tion of "voids" and "Swiss cheese" structures in otherwise 
homogeneous suspensions have been interpreted as mani- 
festations of gas-liquid coexistence |22L l23j . and the small 
lattice spacing of colloidal crystals compared to the one 
expected on the basis of the known density n was inter- 
preted as evidence for gas-crystal coexistence j2J| . These 
possibilities seemed to be confirmed by direct observa- 
tions of (meta-) stable gas-crystal coexistence j^j, and a 
macroscopic gas-liquid meniscus I2fl| although these ob- 
servations were disputed by others |2^.l2^|. 

Despite the ongoing debates due to a lack of experi- 
mental consensus these experimental results, which were 
all performed at low ionic strength with c s in the /iM- 
regime, triggered a lot of theoretical activities to find the 
source of cohesive energy that stabilizes the dense liq- 
uid or crystal phase in coexistence with a much more 
dilute gas phase. The dispersion forces would be the 
first natural candidate to provide the cohesion, but their 
nanometer range is generally considered to be too small 
to dominate over the electrostatic repulsions with a range 



of, 



100 nm at these low salt concentrations. 



It was for instance found that ion-ion correlations, 
which are ignored in the derivation of V2(r), can lead to 
attractive contributions to the pair potential. However, 
the effect is small and too short-ranged for monovalent 
ions at room temperature in water [13|. 

Another avenue of research considered the possibil- 
ity of the breakdown of pairwise additivity due to 
non-negligible effective triplet and higher-order forces. 
Within Poisson-Boltzmann theory the triplet potential 



was calculated, and turned out to be attractive indeed 
thereby suggesting that many-body interactions could be 
the source of cohesive energy. Phase diagrams based on 
repulsive pair interactions and the attractive triplet 
potential indeed showed coexistence of a dilute gas with 
very dense crystal phases (as well as crystal-crystal co- 
existence) [23, E3| j while experimental evidence for the 
breakdown of pairwise additivity was obtained by an in- 
verse Ornstein-Zernike analysis of measured colloidal ra- 
dial distribution functions JjjO], HI H3| > as well as by di- 
rect measurement 0, H3| ■ So although pairwise additiv- 
ity seems to be breaking down at low salinity, it is yet 
questionable whether an approach based on the explicit 
calculation of triplet and higher-order potentials, if feasi- 
ble at all, is very efficient, as convergence is probably slow 
in the regime of strong triplet attractions: there is hardly 
any justification to ignore the four-body potential when 
including the triplet potential changes the phase diagram 
completely compared to the pairwise case. This notion 
was made explicit by a recent simulation study of the 
primitive model (charged colloids and explicit microions) 
that underlies the effective one-component system of 
Ref. [H]: the gas-crystal coexistence that was found 
with included triplet interactions disappeared again in 
the simulations of the multi-component simulation pifil ]. 

An alternative representation of non-pairwise inter- 
actions is based on density-dependent pair-potentials. 
Roughly speaking, this implies that the explicit 
coordinate-dependence of higher-body potentials is 
smeared out to reduce to density dependence in the pair 
interactions. In the case of charged colloids it seems nat- 
ural to modify the form of the screening constant, such 
that not only the background (reservoir) salt concentra- 
tion 2c s but also the finite concentration Zn of the coun- 
terions and the hard-core exclusion from the colloidal vol- 
ume is ta ken into account . For instance, one replaces k 
by k — ^4ttAb(2c s + Zn), J4ttX^2c^_+ ZnVU - rf) or 
similar expressions H El IS I11M El El El El that 
reduce in the dilute limit n — > to n as given by Eq. J5J . 
Often k{n) > k, and one could interpret the resulting 
reduction of the pairwise repulsions due to the more effi- 
cient screening at higher density as an effective attractive 
many-body effect. 

Interestingly, however, a careful analysis of the to- 
tal free energy of the suspension reveals that a density- 
dependent screening constant affects not only the pair- 
interactions but also one-body contributions, such as the 
free energy of each colloid with its "own" diffuse cloud of 

counterions H1E1E1E1E1E1E1EEE3- The Sick- 
ness of this double layer is typically and hence its 
typical (free) energy is of the order of — (Ze) 2 /e(a + /c _1 ), 
i.e. the Coulomb energy of two charges ±Ze at separation 
a + This term lowers progressively with increasing 
n and thus provides cohesive energy, whereas it is an ir- 
relevant constant offset of the free energy if a constant k 
is taken instead of k(n). It was shown that the density- 
dependence of these so-called volume terms could drive a 
gas-liquid spinodal instability at low salt concentrations 



3 



E3- IS 03 j and could hence (qualitatively) explain 
some of the puzzling experimental observations. 



II. HAMILTONIAN, DONNAN ENSEMBLE, 
AND EFFECTIVE HAMILTONIAN 



There are several reasons, however, to revisit the the- 
ories of e.g. Refs.jH H3, E& E3- First of all, they are 
formulated in the canonical ensemble (fixed ion concen- 
trations), which not only obscures its close relationship 
with the classical Donnan theory for colloidal suspen- 
sions (4^.l4fl|. but also unnecessarily complicates the nu- 
merical calculation of phase diagrams as we will argue in 
section II. 

Moreover, and more importantly, the derivation of the 
explicit expressions for the total free energy was perhaps 
not very transparent in Refs. an d may have hin- 

dered extensions of the theory to include, for instance, 
charge renormalization. This nonlinear effect was first 
studied in a cell geom etry j^jj], and, more recently, in a 
jellium-like model[5ll l52l| . In both of those models, the 
nonlinear character of the theory is retained, while its 
complicated multi-centered nature is replaced by a radi- 
ally symmetric structure. The effective colloidal charge 
Z* that appears in the prefactor of the DLVO repulsions, 
is then reduced from its bare value Z due to a tightly 
adsorbed layer of counterions in the vicinity of the col- 
loidal surface. This effect is important when ZXb/o, » 1 
and therefore casts serious doubt|57] 
on the predictions of the gas- liquid and gas-crystal tran- 
sitions in e.g. Refs. pjfl l4"fiL 1471] since large values of Z 
were needed to have the transitions E3| ■ If one now 
realizes that Z* depends on n and K,{n), as was shown in 
e.g. Ref. it is easy to imagine that the volume terms 
are affected non-trivially by charge renormalization simi- 
larly as by the n-dependence of the screening parameter. 
It is therefore important to be able to include this effect 
into volume-term-type theories, and hence to reformulate 
these theories as transparently as possible. 



In order to be able to address all these issues, we re- 
visit here the purely linear screening theory with volume 
terms. Its nonlinear extension to include charge renor- 
malization will be discussed in a forthcoming paper [Hfll ]. 
The present paper is organized as follows. In section 
II we introduce the microscopic Hamiltonian H of the 
colloid-ion mixture and give formal expressions for the 
effective Hamiltonian H for the colloids. In section III we 
calculate H by minimizing the mean-field grand poten- 
tial functional of the ions, whereby explicit expression for 
the density-dependent screening parameter, the Donnan 
potential, and the Donnan effect are obtained as inter- 
mediate results. In section IV we consider the thermo- 
dynamics of the suspension, in particular the free energy 
and the osmotic pressure, with a few interesting cancel- 
ing contributions. In section V we calculate a few phase 
diagrams. We conclude and summarize in section VI. 



We consider a suspension of N identical colloidal spheres 
(radius a, positions Ri, charge —Ze homogeneously dis- 
tributed on the surface) in a continuum solvent of vol- 
ume V characterized by a dielectric constant e at tem- 
perature T. The density of the colloids is denoted 
by n — N/V. In addition there are N + and N_ monova- 
lent point-like cations (+) and anions (— ) present, re- 
spectively, and charge neutrality dictates that N + = 
N_ + ZN. The total interaction Hamiltonian of the sys- 
tem can therefore be written as 



7~L = "ftcc + 'He 



7~L S 



(5) 



where the bare colloid-colloid Hamiltonian Tt cc , the 
colloid-salt Hamiltonian H cs , and the salt-salt Hamilto- 
nian Hss are pairwise sums of hard-core and (unscreened) 
Coulomb potentials. We write 7i cc = J2i<j ^cc(-Sij) with 



/314c (r) 




r < 2a; 
r > 2a, 



and H cs = H c+ +H c - with H c± = £^1 Ef=i K±(|R 



(6) 



r^|), where 



0V c± (r) = < _ZA L! 



r < a; 
r > a, 



(7) 



and is the position of the jth positive (negative) 
micro-ion. The expression for H ss is similar, but without 
the hard-core terms because of the point like nature of 
the ions. 

In principle, the thermodynamic properties of this sys- 
tem could be calculated from the Helmholtz free en- 
ergy of the system F(N, N-, V, T), which is defined by 
cxp(— [3T) = tr c tr+tr_ exp(— (3H). The canonical traces 
are defined as in Eq. <@J. Note that one can ignore the 
explicit N + dependence of T because of the charge neu- 
trality condition. 

Within linearized Poisson-Boltzmann theory and ex- 
ploiting the Gibbs-Bogolvubov inequality, T was explic- 
itly calculated in Refs. ^a, E3- The phase diagram was 
then constructed from T by imposing the usual condi- 
tions of mechanical and diffusive equilibrium, viz. 



(2) „V)\ 



M (n( 1 ),n (1) ) 



= /i(n( 2 \nL 2) ) 
= M-(- (2) ,^ 2) ), 



(«) 



where raW and n_ denote the colloid density and the an- 
ion density in phase i, respectively, and where we intro- 
duced the pressure P = — (dT/dV), the colloidal chemi- 
cal potential \i = (dF/dN), and the anion chemical po- 
tential /!_ = {dT JdNJ). The system (gj of three equa- 
tions for the four unknown densities yielded the phase 
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diagram in the n — n_-plane, for given parameters Z, a, 
and Ab- 

These canonical ensemble calculations were, how- 
ever, numerically rather demanding, since many numer- 
ically expensive evaluations of P(n, /Lt(n, ra_) and 
(i-(n, n_) are needed in the root-finding procedure of 
solving the Eqs. JBJ. For that reason the phase diagram 
of only a few combinations of parameters Z, Ab, and a 
has been studied in some detail. Moreover, the derivation 
of the explicit expressions for T was perhaps not very 
transparent, and may have hindered extensions of the 
theory to include, for instance, nonlinear effects such as 
charge renormalization. And on top of this the (strong) 
connection with the standard description of colloidal sus- 
pensions in terms of a Donnan equilibrium was not made 
in Refs. [Hill. 

It turns out, as we will show in this paper, that at 
least some of these shortcomings and drawbacks of work- 
ing in the canonical ensemble can be lifted by treating 
the anions and cations grand-canonically. For this we as- 
sume the suspension to be in diffusive contact with a di- 
lute reservoir of monovalent anions and cations at chem- 
ical potential p± = fcBTln(c s Aj_), where 2c s is the total 
ion density in the (charge neutral) reservoir, and where 
A± is the thermal wavelength of the cations (+) and an- 
ions (— ), respectively. The colloidal particles cannot en- 
ter the ion reservoir (e.g. because of a semi-permeable 
membrane in an actual experimental setting), and re- 
main treated canonically (fixed TV and V) as before. 
The thermodynamic potential of this ensemble, which we 
will call the "Donnan-ensemble" from now on, is denoted 
F = T — p+N + — fi-N-, and is a function of the vari- 
ables N, V, T, and p±. It is related to the microscopic 
Hamiltonian by the "Donnan partition function" 

exp(-jftF) = tr c Tr+ Tr_ exp(-/3H), (9) 

where H was defined in Eq. © and the grand canonical 
traces are defined as 

Tr± = J2 ex P (/3M±^±)tr± = £ jj-j / dr^ . (10) 

i— 1 i— 1 

Here we have used that exp((3p±)/A± = c s (where the 
factor 1/Aj_ follows from the classical momentum in- 
tegration), and we denoted the microion coordinates 
by r^*. For convenience we will drop the explicit T- 
dependence from now on, and we replace the dependence 
on p± by the reservoir concentration c s . 

Because of the extensive character of F we can write 
F(N, V, c s ) = Vf(n, Cs). The thermodynamic proper- 
ties follow now as p — (dF/dN) = (df/dn) and P = 
— (dF/dV) = np — /, where the derivatives are to be 
taken at fixed c s and T. This implies that the phase- 
coexistence conditions simplify to the two conditions 

P(rc«c 8 ) =/V 2 ),c s ); 



for the two unknown colloid densities at fixed c s , 
i.e. we have prearranged equal chemical potential of the 
ions due to our choice to work in the Donnan ensem- 
ble. This is a considerable reduction of the numerical 
effort compared to the Eqs. JHJ. Note that the me- 
chanical equilibrium condition is equivalent to equal os- 
motic pressure 14 in the two coexisting phases, where 
II(n, Cs) = P(n, c s ) — P(0, Cs) is the suspension's excess 
pressure over the reservoir pressure P(0, c s ) = 2c s k^,T 
(recall that we treated the reservoir as an ideal gas here) . 
This simple relation allows for a rather direct contact 
with experimental measurements of the osmotic pressure, 
as we will also show below . 

Even though our main objective is to calculate 
F(N, V, Cs) as defined in Eq. JjjJ, we will first focus on 
an important and convenient intermediate quantity: the 
effective Hamiltonian H, which depends on the colloid 
configuration {R} and parametrically on the reservoir 
salt concentration c s . It is defined as 

exp(-PH) = Tr+ Tr_ exp(-/3H) 

= exp(-/3H cc )Tr + Tr_ exp(-/3H cs - (3H SS ) (12) 
= exp(-/37Y cc ) exp(-/3fi) 

where, in the last step, we defined the grand partition 
function exp(— (3fl) of the inhomogeneous system of in- 
teracting cations and anions (through H ss ) in the exter- 
nal potential of the colloids (through H cs ). The corre- 
sponding grand potential of this system is f2, which is 
the quantity that we need to calculate in order to find 
the effective Hamiltonian given from Eq. l!12ll as 

H = H CC + fl. (13) 

Once H is known, we can use standard one-component 
techniques to obtain approximate expressions for F, since 
exp(— (3F) = tr c exp(— (3H) is precisely as if F were the 
Helmholtz free energy of a one-component system with 
Hamiltonian H . 



III. THE GRAND POTENTIAL 
A. Density functional 

We will not explicitly calculate from the grand parti- 
tion function of Eq. 11211 . Instead we exploit the frame- 
work of classical density functional theory (DFT), which 
treats an inhomogeneous fluid in an external field at the 
level of the one-body distribution functions (the density 
profiles) (fifllfSlllfi^ : the equilibrium density profiles min- 
imize the (variational) grand potential functional, and 
this minimum is the grand potential. Here we denote 
the density profile of the cations by p+(r), that of the 
anions by p~(r), and the grand-potential functional by 
For notational convenience we do not intro- 
duce a separate notation for the variational and equilib- 
rium profiles, and neither for the grand-potential func- 
tional and its minimum (equilibrium) value. 
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The cations and anions experience external potentials 
U + (r) and U-(r), respectively, due to the Coulomb and 
excluded volume interactions with a fixed configuration 
{R} of colloidal particles. These potentials are explicitly 
given by 



N 



t/±(r)=]>> c± (|R 2 -r|) 



(14) 



where the colloid-ion pair potentials V c ± (r) were defined 
in Eq. (JJ. We can now write the grand-potential func- 
tional within a simple mean-field approximation as 



e 2 f 

Q[p+,p-} = Sl id [/5 + ] + f7 id [p_] + — /drdr 



dr 



(p + (r)U-+(r)+p_(r)t7_(r)), ( 



15) 



where we defined the ion charge density p(r) = p+(v) — 
P- (r) , and where the ideal-gas grand potential functional 
can be written as 

n id [ P± ] = j dr p±(r) (- /i± + k B T (In P± {r)A 3 ± - l] 



k B T y*drp ± (r)(ln 



(16) 



Here we have substituted the identity p± = k B T In c s Aj_. 

The Euler-Lagrange equations 5Q,/5p±(r) — that 
correspond with Eq. 1151 . can be cast, for r outside a 
colloidal hard core, into the form p±(r) = c s exp[=F^(r)]. 
The dimensionless potential </>(r) must then satisfy 
the nonlinear multi-centered Poisson-Boltzmann equa- 
tion jH3 



V 2 0(r) = K 2 sinh0(r) 



ZX 



N 



Sdr-Kil-a), (17) 



where 8(r) is the Dirac-delta. Unfortunately, no analyt- 
ical solution to Eq. l|TYI) is known for the multi-centered 
geometry of interest here. Even solving Eq. 1I17J1 numer- 
ically is far from trivial, and requires a serious computa- 
tional effort 0, El EE El El EH 

For this reason we will first make further approxima- 
tions to the functional, and then perform its minimiza- 
tion afterwards. The main approximation involves the 
expansion, up to quadratic order, of the ideal-gas grand 
potential terms about the, as of yet unknown, ion den- 
sities p±, such that p±(r) — p± are considered to be 
the "small" expansion parameters. This expansion yields 
^id[p±] » O( d [/0±] with 



p ± (\rx^-l)V + \n^ /dr( P± (r)-p ± ) 



P± 



J-Jdr(p ± (r)-p ± ) 



(18) 



In principle, this expansion holds for arbitrary p±, but 
later on we will choose p± to be equal to the average 



ion concentrations in the system, such that /dr (p±(r) — 
P±) = 0) i- e - Vf>± — N± is the number of ions in the sus- 
pension. As will be shown below, this linearization corre- 
sponds to a linearization of Eq. I|17H about 4>(y) — <j> with 
(f> the Donnan potential. This is in line with Ref. [70]. 

It turns out to be convenient, and necessary, to rewrite 
the external potentials U±(r) for the ions such that 
U±(r) = ±V(r) + W(r), where we defined the electro- 
static potential (due to the colloids) V(r) = X^O 1 " — 
Ri|) and the hard-core potential W(r) = J2i w (\ r ~ ^tl)) 
with 



/3v(r) 



and 



0w(r) 



(3v r < a; 

— Z\ B /r r > a, 



(3wq r < a; 
r > a. 



(19a) 



(19b) 



Although we are actually interested in the hard-core limit 
f3v — > oo and f3w n — > oo here, we introduce the (fi- 
nite) hard-core parameters vq and wq here for later con- 
venience. They are necessary and sufficient to ensure, 
within the linearized theory, a vanishing ion density in 
the colloidal hard cores. 

Collecting the results we can write the approximate 
grand-potential functional as 

Q[ P+ , p_] = Q{ d [p + ] + fi>_] + $■ f drdr' 

2e J |r — r'| 

+ /dr (p(r)y(r) + (p + (r) + p.{v))W{r)], (20) 



which is minimized by those (equilibrium) profiles that 
satisfy the Euler-Lagrange equations 

In E± + p± ^ Z£± ± fa) + /3W(r) = 0. (21) 

Cs P± 

Here we introduced the (dimensionless) electrostatic po- 
tential 4>(r), given by 



(r) = A B /dr 



r — r' 



■0V(r). 



(22) 



B. Equilibrium profiles and Donnan equilibrium 

We leave the hard-core parameters vq and wq undeter- 
mined for now, and start the analysis of the Euler- 
Lagrange equations by integrating Eq. (12H over the vol- 
ume. At the same time, we impose that Jdr[p±(r) — 
P±] = 0, i.e. we choose p± such that it is the actual aver- 
age ion density in the suspension. After rearrangement, 
we find that 



p± = c s exp[=F^ - vfiwo], 



(23) 
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where <j> = f dr<fi(r)/V is the spatially averaged electric 
potential, i.e. the Donnan potential. Since global charge 
neutrality imposes that p+ — p- = Zn, we can conclude 
from Eq. 12.H that the Donnan potential satisfies 



sinh ( 



Zn 

-7T~ ex P [VPWo \, 
2c s 



(24) 



which reduces to the usual Donnan expressions in the 
point-colloid limit 77 — » [illlll. Combining Eq. d 
with yields 



P± 



\ [yJ{Znf + (2c s ) 2 exp(-277/3 W o) ± Zn) , (25) 



which explicitly relates the salt concentration in the sus- 
pension to the colloid density and the salt reservoir con- 
centration, provided the parameter wq is known. 

Using these relations for </> and p± we consider a specific 
linear combination of the Euler-Lagrange equations, and 
rewrite Eqs. (l2ll) as 



2 = - 2 (PW(r) - n[3w Q ) ; (26a) 

P = -{p+ + P-){<i>{*)-4>) 



p+0) , p-( r ) 
p+ p- 

P(r) 



p((3W(r)-r)0Wo), 



(26b) 



where we defined the short-hand notation p — p+ — p- — 
Zn for the overall ionic charge density. This particular 
linear combination was chosen, because (i) the charge 
density is the physical quantity of interest here, and 
(ii) the electric potential is decoupled from the "charge- 
neutral" density. 

It is straightforward to solve the "hard-core" linear 
combination, Eq. (|26a|) . Imposing that p + (r)/p + + 
P-(r)/p~- = within the hard-core of any of the col- 
loids (i.e. wherever W(r) = wq) yields a value for the 
hard-core parameter, 



(3w 



1 



1-r,' 



(27) 



whereas outside any of the colloidal hard core positions 
we have 



P+(r) , P-(r) 



Pa 



P- 



1 



>1 



(28) 



The solution of the "charge" linear combination, 
Eq. l|26hjl . is most straightforwardly found by Fourier 
transformation. For an arbitrary function /(r) we 
define and denote the Fourier transform as /k = 
Jdr /(r) cxp(ik • r). One easily checks from equa- 
tion lt26hTl that 



(29) 



p k =(2tt) 3 {p + (p+ + p-)4> + pvpwo} J(k) 
- (P+ + P-)<k* - pW k , 
where we have from Eq. l|19b|l that 

Airawo ( sin(fea) 



Wu 



k 2 



ka 



\ N 

cos(fca)) ^V k - R ^, (30) 
' 3=1 



and from Eqs. |2J an( ^ P^ajl that 

A N 

Pk 47ra v 
<^k =4ttA b - -T2- 2^ ex P(* k • R j) 

\( a , v Xb \ 1 a sinfca 
x < pv + Z — cos/ea — pvo — ; 

\ a J ka 



(31) 



Equations H29J1 and 13 H are two linear equations in the 
unknowns 4>k an d Pk, which can be solved straightfor- 
wardly. Fixing the remaining hard-core parameter vq to 



a y K ^B ,„ P+~ P- 

(3v = -Z _ + /3w - — — - 
1 + Ka p + + p- 

we find that the charge density is given by 
Pk =(2tt) 3 j 



(32) 



+ {p+ + P-)(t) 



r 5(k) 



+ 



1 — n v T ' T J k 2 + k 2 
Z cos ka + f sin ka — -- (3 ' ] ' 



1 + Ra 1 + k 2 j~k 2 



where the effective Debye screening parameter is defined 



K =^47rA B (/3+ + p-) 

= ^/4^i/(Zn) 2 + (2c s exp[-77/(l - n)}) 2 



(34) 



Here we used Eqs. ll25|l . and ll27l in rewriting 

the first into the second line. Note that the factor 
exp[— 77/(1 — n)] that appears in Eq. 13411 can be accu- 
rately represented by (1 — 77), with a relative deviation 
less than 0.01 for 77 < 0.1 and less than 0.1 for 77 < 0.35. 

The first term in expression 13311 is of the form oc 
/c 2 <5(k) and does not contribute to the charge density Q35H ■ 
However, we will see below that this term does in fact 
contribute to the grand potential, as this involves the 
Coulomb energy oc Jdkp-^/k 2 . 

The real space representation of the charge density 
is a multi-centered sum p(r) — J2iPi(\ r ~ R-il), where 
the one-particle density profiles (the "orbitals") have the 
usual DLVO formal H(|: 







Pi( r ) — { Zk 2 exp(«a) exp(— nr) 



4vr 1 



r < a; 



r > a. 



(35) 



We note that the vanishing of pi(r) inside the colloidal 
hard core is a direct consequence of our particular choice 
for vq given by Eq. 113211 : other choices for vq would have 
yielded a finite ion charge density inside the hard core. 
Note also that the multi-centered charge density p(r) is 
not vanishing within the hard cores, since the exponential 
tail of the orbital centered around colloid i penetrates the 
hard core of all the other colloids j ^i. 

By inserting Eq. il2"Tjl into il25ll , explicit expressions for 
the average concentrations p± of ions in the suspension 
are obtained as a function of the colloid density n, colloid 



7 




(b) 



Figure 1: Total ion concentration p+ + p- (a) and concentration of added salt p_ (b) as a function of the colloid packing 
fraction for different reservoir concentrations, using the expressions of Eq. 125|l and l|27jl . The colloidal charge and radius are 
Z — 50 and a = 21. tarn, respectively, and the solvent is ethanol at room temperature such that Ab = 2.37nm. This matches 
the parameters from the experiments by Rasa et al [71 ]. 



charge Z, and the reservoir concentration 2c s — this was 
already used to obtain Eq. p3|). These expressions re- 
duce, in the limit of point-like colloids (for which 77 = 0) 
to the standard expressions for the Donnan effect |48l 

This effect is illustrated in Fig. where we plot the 
total ion concentration p+ + p_ in (a) , and the concen- 
tration of added salt 2p- — p + + p- — Zn in (b), on 
the basis of our expressions for p± . The param eters are 
close to those of the experiments hv lRasa et alJ [z3,[z3|: 
Z = 50, Ab = 2.37nm, and a = 21.9nm. The reser- 
voir salt concentration equals the 77 = limit of each of 
the curves, and the crossover from the low-77 plateau to 
the high-77 linear part corresponds to the crossover from 
added-salt dominance to counterion dominance. Note 
the expulsion of added salt back into the reservoir at 
high 77 in (b) . An important aspect of these intermedi- 
ate results is that the scre ening parameter R increases 
with n essentially oc V Zn in the counterion-dominated 
regime (which may occur at packing fractions as low as 
77 ~ 1(T 4 if c s ~ 3/uM). 

As we have now solved the Euler-Lagrange equations 
(12H for the two linear combinations p+(r)/p + + j o_(r)/ / o„ 
and /?+ (r) — (r) , it is straightforward to disentangle the 
equilibrium profiles and obtain the profiles p± (r) of the 
two ionic species separately. 

It is important to realize, however, that these results 
depend on the particular choice that we have made for the 
hard-core potentials in Eqs. Hlflajl and ijlflhjl . Different 
choices for these hard-core potentials lead to other, non- 
equivalent minima of the grand-potential. For instance, 
instead of U±(r) = ±V(r) + W(r), we could have consid- 
ered the choice U±(r) = ±V(r) + 2p T W(r)/(p + + p_), 
which, with (3vq = —Zk\b/(1 + Ra) and /3wq = 1/(1 — 77) 



would lead to a vanishing p\{r) and p + (r ) / p + + p^ (r ) /p_ 
inside the hard cores. This choice was actually made in 
Refs. miE3; an d leads to similar, but not identical re- 
sults — see Appendix IbI 



C. The minimum of the functional 

In Appendix El we derive the equilibrium grand poten- 
tial fl by insertion of our solution of the Euler-Lagrange 
equations into the functional. The effective interaction 
Hamiltonian H — Ti cc + then takes the form 



H({R},N,V,c s 



N 
i<j 



;n,c B ). (36) 



The first term is independent of the colloidal co- 
ordinates Ri, and is called the "volume term" as it is 
a density-dependent, extensive thermodynamic quantity 
that scales with the volume of the system. The sec- 
ond term of Eq. l-'Uit is a pairwise sum that does de- 
pend on the colloidal coordinates (and on the density n). 
For later convenience we decompose the volume term as 
$ = $d + $0; with the so-called "Donnan" term defined 
by 



/3£d 
V 



E 



p± 



ln^-1 



and the other term by 
1 (Zn) 2 



V 



2p+ 



r] 2p + p- n Z 2 k\b 



1 



2 1 



(37) 



(38) 
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Tl=0.01 

ti=0.001 
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Figure 2: The factor A of Eq. I|40fl as a function of Ra and 
r\. For all parameters, A <C 1, so it can safely be neglected in 
Eq. (O. 



U(r) = V(r) + 2p T W(r)/(p + instead of the defini- 

tion used here. This latter choice does not affect any of 
the other volume terms, but does involve another choice 
for vq and wq, and does change the expression of <f> (see 
Appendix El . For these two reasons we set A = in the 
remainder of the paper. 

The so-called volume terms 3>d and $ are very sim- 
ilar to their canonical counterparts that were derived in 
Ref. ^jj- The main difference is that the present vol- 
ume term includes the term — Jdr (/x+p+(r) +/i_p_(r)) 
due to the grand-canonical character of our calculations. 
This leads to another difference, since one should now 
view the Hamiltonian H as a function of n and the reser- 
voir salt concentration c s , i.e. one should take the de- 
pendence of <f> on p± and R as a dependence on c s and n 
through Eqs. I.'t4ll and 12511 . It is the nontrivial (and non- 
linear) dependence of /3^u/V and (3<&o/V on the colloid 
density n, at fixed c s , that is responsible for interesting 
thermodynamic effects, as we will see below. 



In Section TlV Bl below, we will see that $d, which takes 
the form of ideal-gas contributions, accounts for the Don- 
nan equation of state (except for the colloidal ideal gas 
contribution); hence the nomenclature. The term $ ap- 
pears as an electrostatic (and hard-core) free energy con- 
tribution. This separation is slightly misleading, how- 
ever, since the two terms both depend on n and Z 
through the expressions ll25)l and iJHU, which stem from 
the Donnan potential II24H and hence from the balance 
between electrostatics and entropy. 

The effective pair potential between the colloids, 
V(Rij), that appears in the second term of Eq. J3EJ), is 
given by 



r < 2a; 
r > 2a, 



with the DLVO-charge given by Z > = Z exp(na) / (1 
Ra), and with the parameter A defined by 



A = inpwo— ((1 + Ra) 2 e' 2Ra + (Ra) 2 - l) 



(40) 



The effective pair interaction V(r) is very similar to the 
traditional DLVO potential V 2 (r) of Eq. QJ, but with 
two important differences. The first difference involves 
the screening parameter R in V(r), which is to be con- 
trasted with the reservoir screening parameter n in V2 (r) . 
The second difference is that the amplitude of V(r) is en- 
hanced compared to Vzir) by a factor (1 + A). This can 
be traced back to our particular choice of linear com- 
binations of density profiles that we used to solve the 
Euler-Lagrange equations. 

In Figure we plot A as a function of the screening 
parameter Ra for several values of the packing fraction 
77. The plot shows that A <C 1 for essentially all pack- 
ing fractions of interest here. Moreover, one can also 
show that A = if the hard-potentials are defined as 



IV. THERMODYNAMICS 

A. Free energy 

As we have now found the functional form ll.'ifill for the ef- 
fective Hamiltonian of the colloids, we are ready to calcu- 
late the corresponding free energy F(N, V, T, c s ) defined 
just below Eq. ll.'tll . From this the other thermodynamic 
quantities follow. Since the volume terms in l-'Kit are 
independent of the coordinates of the colloids, we can 
factor out their Boltzmann weights and write 



exp(— (3F) = cxp[— /3<t>]tr c exp 
This can be rewritten as 



N 



-pJ2V(Ri, 



F = $ D + $0 + F id + F e: 



(41) 



(42) 



with <3?d and <f>o defined in Eq. Q,'i8[l . with the colloidal 
ideal-gas free energy 



F id = Nk B T(\n(nV) - l), 



(43) 



and where F exc is the non-ideal (excess) free energy due 
to the colloid-colloid pair interactions (l3Qf) . Here we cal- 
culate Fpjr yariationally, using the Gibbs-Bogoliubov in- 



equality [74L I75L 1 7(1 177| . This inequality states that the 
excess free a so-called reference system 

of volume V that contains N particles with any pair in- 
teraction V^(Rij), satisfies 

^ X c<^s f) + (E( y (^)- F(ref) (^))) < ( 44 ) 

\i<3 I ref 

where (. . . ) re f denotes a thermodynamic average that is 
to be evaluated in the reference system. The key idea is 
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to use a reference pair potential with a free parameter 
with respect to which the right hand side of Eq. lliHl can 
be minimized; the minimum is then the optimal estimate 
for the free energy F exc of interest. We use two different 
reference system to calculate the free energy of fluid and 
crystal phases, respectively. 

For the fluid phase we use a hard-sphere reference sys- 
tem, with the hard-sphere diameter d as variational pa- 
rameter. Introducing the effective hard-sphere packing 
fraction £ = (ir/6)nd 3 , we can write 



fexc . f 4£ - 3£ 2 

= mm< -rr 

Nk B T d \(l-0 2 

n , 
+ -4* 



drr 2 g d (r;OPV{r) 



(45) 



where the first term is Carnahan-Starling expression for 
the hard-sphere free energy [74ll78T |. and where gd(r; £) is 
the radial distribution function of a fluid of hard spheres 
with diameter d and packing fraction £. We approximate 
g d (r;£) by the Verlet-Weis corrected Percus-Yevick ex- 
pressions [lillzll- This allows for an analytic evaluation 
of the integral in Eq. "J"}, since the Yukawa form of V(r) 
turns this integral into a Laplace transform of rgd(r;£), 
for which accurate expressions have been derived on the 
basis of Pade fits in Refs. (sfT . Isij . The minimization 
with respect to d is then numerically performed straight- 
forwardly. Note that such a minimum indeed exists, as 
the excess free energy of the hard-sphere reference sys- 
tem becomes infinitely large in the limit of large particle 
sizes. Because the particles in our actual system have 
also a hard core of radius a, we impose that d > 2a. 

As a reference system for the solid phase, we use 
N classical Einstein oscillators [82L Isflj on an FCC lat- 
tice The Einstein frequency loe plays the role of the 
variational parameter used to minimize the right hand 
side of Eq. 14411 . For the thermodynamic average of the 
Yukawa intera ctions in th is system, we use the expres- 
sions found by IShih et all [84j . 

So far we have only considered FCC configurations for 
the solid phase, but there is no principal problem to gen- 
eralize this to other structures such as BCC j^3. Con- 
sequently, we only consider gas-liquid, fluid-FCC and 
FCC-FCC phase equilibriums in this paper. 



B. Osmotic pressure 

The osmotic pressure IT = P — 2c s kT of the suspen- 
sion under consideration follows from P = —dF/dV at 
fixed N and c s . We can therefore use our expression for 
F given in Eq. *32|) to obtain II(n, c s ) explicitly as 



n = n D + n + n id + n e 



(46) 



with n D = -2c s fe B T- {d&o/dV), the Van 't Hoff (ideal- 
gas) contribution n id = — (dF id /dV) = nfc B T, the ex- 
cess pressure II exc = — (dF exc /dV), and the remaining 



term TIq = — (d^o/dV). Explicit general expressions 
for LTd, LTo, and LT exc can of course be given on the ba- 
sis of Eqs. (l38l) . and e.g. 1(4*51) . respectively, but it 
turns out to be instructive to focus on these expressions 
in the limit of point-colloids with radius a — (such that 
T] = 0): this reduces the algebra and allows for an in- 
teresting illustration of cancellations of some of the elec- 
trostatic contributions to the osmotic pressure LT. We 
stress, however, that we used the full expressions in our 
numerical calculations presented below. 
In the point-colloid limit we have 



/3II D = - 2c s + p+ + p- 

= - 2c s + y/(Zn)* + (2c) 2 

(Zn) 2 



(47) 



0(n 4 ) 



Zn -C 2c s 



4c s 

Zn - 2c s + 0{c 2 s ) Zn > 2c s , 



and a little tedious but straightforward algebra yields 

{Zn) 2 



(3Uo 



4c. 



0{n 6 



Zn <C 2c s 



-bn 3/2 + (D(c 2 ) Zn-»2c 



(48) 



with a coefficient b = V 1 Ztx\bZ 2 \b /2. We focus first on 
the low-density /high-salt regime Zn <g; 2c s , and then on 
the opposite regime. 

The expressions 114 711 and 1481 show a cancellation of 
the dominant term in the regime Zn -C 2c s , such that in 
this regime II ~ Ilid + II exc , i.e. the pressure is actually 
the pressure of the effective one-component system de- 
scribed by the pairwise screened-Coulomb Hamiltonian. 
Interestingly, however, one can also write the virial ex- 
pansion (3F exc /V = B2{R)n 2 + Q(n 3 ) in this regime, 
where the second virial coefficient 



B 2 (K) 



exp[ 



-pV(r)]), 



(49) 



with the colloidal pair potential V(r) defined in Eq. 13911 . 
In the limit of weak interactions, the exponent in Eq. Q49|l 
can be linearized with the result that B 2 = Z 2 /4:C S for 



point-colloids. This means that /?II e 



(Zn) 2 /Ac 



/3IId, and hence that the pressure can also be approx- 
imated by the Donnan expression IT ~ ITid + LTd- In 
other words, on the basis of this simple analysis one ex- 
pects "reliable" results for the pressure (and hence for the 
thermodynamics) in the regime Zn <C 2c s by taking ei- 
ther the full four-term expression 14611 for LT, or the two 
two-term expressions Ilid + LTd and Ilid + n exc , but not 
any other combination. This will be confirmed by our 
numerical results below. 

The situation is bit more complicated in the opposite 
low-salt regime 2c s <C Zn, since then (i) no cancella- 
tions take place and (ii) the virial expansion for F exc 
breaks down because of the long-range character of the 
unscreened-Coulomb interactions. 
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Figure 3: Eq uation of state compared to the experiments of 
Ra§a et al|72|. described by the parameters Z — 32 for the col- 
loidal charge, a/A B = 96.7 for the radius-to-Bjerrum length 
ratio, and 2c s = 16 fjM for the reservoir salt concentration. 
Shown are the experimental data (crosses), a one-component 
DLVO system (dotted), the Donnan theory (dashed), and 
our full linear theory (solid line). Two of the three theo- 
retical curves describe the experimental curves accurately for 
rj < 0.003, the Donnan theory is less accurate although still 
qualitatively reasonable in this regime. 

As a simple approximation for highly charged particles 
(specifically, particles for which Z 2 \^/a 3> 1), the pair 
correlation function gd(r) can be set to gd(r) = 1 for 
r > in -1 / 3 and to otherwise. One can then show 
that the lowest order contribution to the excess pressure 
takes the form /3II CXC = -6'n 4 / 3 with b' = ttZ 2 \ b /12. As 
a consequence we find the asymptotic low-salt result 

m = (1 + Z)n - b'n 4/3 - bn 3/2 , (50) 

which contains Donnan, colloidal pair, and Debye- 
Hiickel-like contributions. The prefactors of the frac- 
tional powers would change if a proper Giintelberg charg- 
ing process would have been performed but the 
present analysis is good enough to capture the spinodal 
instability that is now well-known to be realistic for prim- 
itive model systems at sufficiently strong coup ling (low 
enough temperature) jH El Hi El S HI On 
this basis one could expect that the present theory pre- 
dicts phase-separation in low-salt colloidal suspensions. 
Within the full theory for F we indeed find this phe- 
nomenon in the next section. 

We now illustrate our results for the osmotic pressure 
by numerically comparing the theoretically predicted val- 
ues to experimental measurements in Figure The ex- 
perimental system is an ethanol suspension of colloidal 
silica-spheres, for which H(n) was determined by inte- 
gration of the measured density profile in sedimentation 
equilibrium [73 , l73j . The system parameters are Z = 32, 
2c s = 16 /xM, Ab = 2.38 nm and a = 21.9 nm. Since 
ZX B I a w 3 we do not expect too much charge renormal- 
ization, and as Zn/2c s Rj 0.7 at the highest density con- 



sidered here (rj = 0.01), this experiment is expected to be 
in the high-salt regime where not only the full expression 
ijlfil) for II but also both the one-component expressions 
11 ~ IIid + n exc and the Donnan expression II ~ Ilid+IlD 
are expected to "work" with reasonable accuracy. 

This is to some extent confirmed by Fig. where the 
measured osmotic pressure is in quantitative agreement 
with two of the three theoretical versions at low pack- 
ing fractions r\ < 0.003 or so; the Donnan pressure is 
less accurate. At higher densities the different theoreti- 
cal curves deviate from each other (and from the exper- 
iment), with the one-component result Ilia + II exc be- 
ing closest to the actual experiment. A word of cau- 
tion is appropriate here, however, since recent work by 
Biesheuvel indicates that charge regularization is relevant 
in the present system, i.e. the bare colloidal charge Z is 
not a constant but decreases with density, where signifi- 
cant deviations of the low-density charge is predicted for 
rj > 0.002 [93|- This is rather precisely the regime where 
the theories begin to deviate from the experiment. This 
issue will also be addressed in more detail in future work. 

From the fact that the one-component osmotic pres- 
sure II = Ilid + II exc describes the experimental data 
rather accurately, one may conclude that the experimen- 
tally found "inflated" profiles of Ref. 0| need not neces- 
sarily be described by theories such as those of Ref. 
where a three-component mixture (cations, anions, and 
colloids) in gravity gives rise to an ion-entropy-induced 
self-consistent electric field that lifts the colloids to higher 
altitudes than expected on the basis of their mass. The 
equation of state suggests that an alternative description 
could be given, based on hydrostatic equilibrium of an 
effective one-component system of colloidal spheres with 
pairwise screened-Coulomb repulsions only. The latter 
picture is not in contradiction with the existence of the 
electric field, since the density variation with height im- 
plies a variation of the Donnan potential with height 
through Eq. pljl. The two pictures are, in this sense, 
merely two sides of the same coin, at least on length scales 
beyond which the local density approximation applies 
that underlies the one-component theory. On smaller- 
length scales the source of this electric field involves de- 
viations from local charge neutrality, which cannot be 
explained by hydrostatic equilibrium and a bulk equa- 
tion of state alone. 

The other system for which we calculate the osmotic 
pressure is one of the systems that Linse studied by 
Monte Carlo simulations in Ref. [95] . This system is free 
of added salt, contains colloids with a charge Z = 40 and 
a radius-to-Bjerrum length parameter of a/As = 22.5 for 
monovalent ions (in the notation of Ref. j^al the cou- 
pling parameter is Tu — 0.0445). The simulated re- 
sults are shown in Figure 0] together with several ver- 
sions of the present theory. It is clear that the major 
contribution to the osmotic compressibility factor orig- 
inates from the pressure n id + II D Rj [Z + l)n, which 
exceeds the one-component combination IL^ + ll exc by 
at least an order of magnitude. The decrease of 0IL/n 
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Figure 4: Equation of state (compressibility factor) com- 
pared to the computer simulation of Ref. where Z = 40 
and a/Aa = 22.5. Shown are the simulation data from [flfil ] 
(crosses), the pressure II from our full linear theory (solid 
line), the approximate low-salt expression IfiUI (dotted), and 
the pairwise one-component result Ilid + n ex c (dashed). The 
theoretical curves are based on the reservoir salt concentration 
c s = 10 -15 M, which is low enough to ensure an essentially 
vanishing coion concentration in all state points shown here. 



Figure 5: Phase diagram for a colloidal suspension as a func- 
tion of colloidal packing fraction 77 and reservoir salt concen- 
tration c s . The colloidal radius and charge are a — 326 nm 
and Z = 7300, and the solvent is water at room temperature 
such that the Bjerrum length is Ab = 0.72 nm. The solid 
lines denote fluid-solid and solid-solid binodals, and the dot- 
ted line shows the underlying metastable gas-liquid binodal. 
The fluid-solid-solid triple point is denoted by x, and the 
solid-solid critical point by □. 



f° r V ^ 0-02 is due to the contribution Ho. Our calcu- 
lated pressure describes the simulation data quite well, 
showing that volume terms may have a pronounced ef- 
fect on the thermodynamic properties of low-salt suspen- 
sions, while the pairwise DLVO-picture without volume 
terms breaks down qualitatively. We note, finally, that 
the limiting expression 15(H for the pressure in the limit 
for point-colloids can be seen to catch the low-density 
negative curvature of (311/ n with n as predicted by the 
full theory and the simulations, but not the increased 
stability at higher n. 



V. PHASE DIAGRAMS 

From the free energy per unit volume f(n) = F/V at 
fixed c s , we calculate the chemical potential and the pres- 
sure, and we impose the usual condition of thermody- 
namic equilibrium 11 111 to find a phase-equilibrium. We 
already mentioned that this is numerically much less in- 
volved than in the canonical calculations of e.g. Ref. |47j . 
where the set (jHJ is to be solved. We merely illustrate the 
feasibility of these calculations here by showing two phase 
diagrams, for a particular Z, a, and Ab- In forthcom- 
ing publications we will fully exploit the relative simplic- 
ity of the grand-canonical formulation of the theory by 
"scanning" the full parameter space (Z, As/a), including 
a generalization of the present theory to include charge 
renormalization jo?! ]. 

The first set of parameters that we consider is Z = 
7300, a — 326 nm, and Ab = 0.72 nm, which corresponds 
to the experiments of Ref. j2^| . The phase diagram that 



follows from the present theory is displayed in Figure 
and shows phase coexistence with a considerable or large 
density gap at c s < 20 /zM, and only a very small den- 
sity gap at higher c s . At a salt concentrations of about 
23.8 /LtM, a liquid-solid-solid triple point occurs (denoted 
by x in Fig.OjJ, and at 26.2 jjM. a solid-solid critical point 
is located (denoted by □ in the figure). Although some- 
what difficult to see in this picture, there is no lower 
critical point. 

The phase diagram of Fig.0is pretty similar to the one 
calculated in Ref. using the canonical version of the 
theory [47j |. but with a few substantial differences. The 
canonical theory, for instance, does not find any solid- 
solid coexistence, nor does it find a triple point for these 
parameters. Also the canonical theory predicts a lower 
critical point, while the current grand-canonical version 
of the theory does not. Despite this differences the main 
phenomenon is shared that at low salinity c s < 20 /jM a 
density gap opens up. 

The physical mechanism for this demixing transitions 
into a dilute and dense phase is identical to what was 
explained in Refs. 0,1221 the self energy of the dou- 
ble layers as represented by the third term of the volume 
term $0 in Eq. I|38|) drives a spinodal instability at low 
enough c s , even though the pair interactions are purely 
repulsive. The underlying physical mechanism is that 
the cohesive energy that stabilizes the dense phase stems 
from the compression of the double layers thickness k" 1 
upon increasing the colloid density: this effect brings the 
charge in the diffuse double layer closer to the oppositely 
charged colloidal surface. This mechanism is very sim- 
ilar to the one that causes gas-liquid demixing in the 
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il 

Figure 6: Phase diagrams for the pa rameters of the exper- 
iment hv iMonovonkas and Oastl |lfl(1 |. The Bjerrum length 
for this system is Ab = 0.72 nm and the partical radius 
is a = 66.7 nm. The data points plotted here correspond to 
the sam ples for which a fluid-solid coexistence was observed 
in Ref. jlflflj . The solid line is a phase diagram for a one- 
component DLVO system, the dashed line denotes the phase 
boundaries for our full linear theory for a charge of Z = 1200. 
The dotted line gives the phase boundary for the same full 
linear theory, but with a lower charge Z = 720. 



restrictive primitive model according to Debye-Hiickel 
theory jUii. 

A word of caution is appropriate here: given that 
ivAe/a — 16, one expects a substantial renormalization 
of the charge within nonlinear Poisson-Boltzmann theory 
for this system, and hence a reduction of the tendency to 
demix. Whether or not this mechanism for phase separa- 
tion remains strong enough to yield a big density gap in 
the phase diagram if charge renormalization is taken into 
account, will be investigated in a future publication [5fll |. 

The second phase diagram that we present 
her e is for the parameters of the experiments 
of iMonovoukas and Gastl [l0rT |. where Z = 1200, 
a = 66.7 nm, and Ab = 0.72 nm. These parameters 
were chosen because the experiments reveal a significant 
density gap, by a factor of three, between the coexisting 
fluid and solid phases at salinity of the order of 10 /LtM, 
Such a large density gap cannot be explained by the 
DLVO pair potential alone, and hence we investigate 
here to what extent the volume terms may account for 
this effect. 

The phase diagram, shown in the rj-p- representation 
in Fig. El shows the experimental points and three fluid- 
solid binoda ls ba sed on the present theory. As the re- 
sults of Ref. |l00l | seem to be independent of the concen- 
tration of added salt for salt concentration lower than 
approximately 8 jjM, we assumed an extra background 
salt concentration of 8 /iM for the experimental points. 
Note that this representation of the phase diagram, with 
the vertical axis representing the concentration of added 
salt instead of the reservoir concentration, is such that 



the tie-lines (which have been omitted for clarity) are no 
longer horizontal as in the r]-c s representation, but in- 
stead tilted to lower p_ at higher i] due to the Donnan 
effect (see also Fig.^l. 

The first binodal in Fig. UJ is the one based on the 
ideal and excess part of F only, i.e. we assume that $o 
and <I>d vanish (or more accurately: the volume terms 
are assumed to be merely linear in N and V and do 
therefore not affect the phase diagram). Although this 
binodal gives a fair representation of the experimental 
points (probably this is how Z — 1200 was chosen), they 
do not capture the large density gap. The second binodal 
is based on the full expression for F, including the vol- 
ume terms, with Z = 1200. We find an enormous density 
gap that is much larger than experimentally observed, 
and that extends to unreasonably high salt concentra- 
tions. The third binodal is also based on our expression 
for F with volume terms, but now for a smaller charge 
Z = 720. Interestingly, this choice gives a density-gap 
at fluid-solid coexistence in the right salt concentration 
regime, but the magnitude of the gap is yet much big- 
ger than experimentally observed. The reduction of the 
charge from Z = 1200 to Z = 720 may give a rough idea 
of the effect of charge renormalization, and shows that 
this nonlinear effect reduces the tendency to demix con- 
siderably. Theories for charge renormalization |5fl l5^ | 
show that the renormalized charge is actually not a con- 
stant but depends on the screening parameter and the 
density; the present value Z = 720 corresponds to the 
dilute limit value at kcl w 0.8, i.e. at c s 10 /LtM, and 
is fixed here for simplicity. We expect this to be a rea- 
sonable lower limit for the renormalized charge in the 
region in which the phase-separation occurs. Also this 
system will be investigated within a nonlinear version of 
Poisson-Boltzmann theory in a future publication [59] . 



VI. CONCLUSION 

We have reformulated and re-derived the volume term 
theory for suspensions of charged colloids HHI3. Our 
present derivation should be more transparent than the 
original one, for instance because we can now avoid the 
extra parameter A that regulates the Coulomb potential 
from 1/r to exp(— Xr)/r with A — > only at the end of 
the calculation. Moreover, the presently derived expres- 
sions should be easier to use in numerical calculations of 
thermodynamic properties and phase diagrams, because 
the ions are treated grand-canonically instead of canon- 
ically, thereby assuring equal chemical potential of the 
ions from the outset. Moreover, a direct connection with 
Donnan theory is now made, with explicit expressions 
for the Donnan potential and the ion concentration in 
the system. In future publications we will fully exploit 
the computational advantages and extend the theory to 
include charge renormalization. 

We derived analytic expressions for the osmotic pres- 
sure in the point-colloid limit for both the low-salinity 
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and high-salinity limits. The low-salinity limit of the 
pressure was shown to correspond to the Donnan ex- 
pression, while in the limit of high salt concentrations 
the traditional DLVO results are recovered. The present 
full theory interpolates between these results, and gives 
a good account of measured and simulated osmotic pres- 
sures in both regimes. 

We also calculated two phase diagrams. The first one 
matches the parameters of Ref. and shows a similarly 
large phase-instability at low salinity, although there are 
also a few substantial differences. The second phase di- 
agram matches the pa rameter s of the experiments by 
Mon ovoukas and Gast of Ref. 100], where an anoma- 
lously large density gap at fluid-solid coexistence was re- 
ported. Interestingly, the present theory does predicts a 
density gap at fluid-solid coexistence, but its magnitude 
is much larger than experimentally observed. We stress, 
however, that these phase diagrams are calculated in a 
regime where charge renormalization cannot be ignored. 
The relative transparency of the present derivation al- 
lows to systematically include this nonlinear effect into 
the theory, as will be shown in a future publication [Hfl| . 

The linear theory described in this paper already 
shows, however, that volume terms can affect the os- 
motic pressure of low-salt suspensions qualitatively, also 
in regimes where charge renormalization and other non- 
linear effects are not expected to be important. 
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where we inserted the Fourier transform p k of p(r) 
from l(33jl . and the Fourier transform Vk of V(r), which 
is given by 
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Appendix A: THE GRAND POTENTIAL 



In this Appendix, we derive the equilibrium grand po- 
tential Q. We show that upon insertion of this grand 
potential into Eq. Ill -ill , the effective Hamiltonian can be 
cast into the form specified by Eqs. H3fi|l -H39 |l . 

In the framework of Density Functional Theory, the 
equilibrium grand potential is given by the minimum of 
the functional Vl[p + ,pJ\ of Eq. H20I1 . This minimum is 
found by inserting the Euler-Lagrange equation l|2"T|l into 
the functional. This leads to the following expression for 



The factor A/2 in the fourth term on the right hand 
side of Eq. l|A2ll is caused by the expulsion of microionic 
charges from the colloid cores, and is given by Eq. lHf))l . 
Note that the first and second term of Eq. l|A2ll result 
from the oc /c 2 <5(k) term in equation 13311 that did not 
contribute to the charge density. 

In a similar way, the "hard-core" part of the grand po- 
tential (I A If) is evaluated as 
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where the Fourier transform of W(r) is given by Eq. l(3"0"jl 
and where we used that W (r ) [p + (r) / p + + p_ (r) /p_ ] = . 

Substitution of Eqs. l|A2ll and l|A4ll into the grand po- 
tential l|Alj) leads to 
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where the "volume term" $ = $d + $0 is given by 
Eqs. {23 and {HHt - 

Gathering Eqs. {J3J, @ and l|A5|) . we find that the 
effective Hamiltonian can be cast into the form given by 
Eqs. Ip|-{39|. 



Appendix B: ALTERNATIVE HARD-CORE 
TERMS 



We have already mentioned that, in this paper, we use 
a slightly different definition of the hard-core parameters 
/3i>o and 0wn than were used by Van Roij and Hansen 
in Refs. |43,|43|. In this Appendix we make this state- 
ment explicit and calculate, within the grand-canonical 
scheme of this work, the effective Hamiltonian us ing the 
definition of the hard-core potentials of Refs. |4fitl47l|. 

In contrast to the presently used definition U±(r) = 
±V(r) + W(r) for the micro-ion-colloid interactions, as 
outlined just above Eq. I|19all . Van Roij and Hansen used 
the definition 



U±(r) = ±V(r) 



2p=, 
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where the potentials V(r) and W(r) are defined in 
Eqs. (I19a,|) and l|19hj) . With this definition, the grand 
potential becomes 
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W{v), 



where the ideal-gas functionals /3f2id[/0±] are defined in 
Eq. ®, 

The corresponding Euler-Lagrange equations are then 
given by 



ln P ± + P ± ir 1 ^ ±Hr) + 2j ± pW(r 1 

Cs P± P+ + P- 
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By integrating these equations over the system volume, 
and using the condition for global charge neutrality, we 
find that the average densities p± are identical to those 
given in Eq. 12oll . The Donnan potential <j>, however, is 
not given by Eq. <24ll anymore, but by 



sinh 



Zn 
2qT 



n(3w - 



Zn 



instead. Although this expression also reduces to the 
usual Donnan expression in the limit n — ► 0, it is physi- 
cally less satisfactory than the result we found in Eq. ll24*jl , 
as at high 77 its sign can become different from that of 
the colloidal charge — Ze. 



To calculate the density profiles, we take the fol- 
lowing linear combination of the Euler-Lagrange equa- 
tions iB? 



=2 (1 - /W(r) + vPwq) ; (B5a) 
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Note that, due to the different definition of U±(r), the 
hard-core potential W(r) is now totally decoupled from 
the charge density p(r). 

Eq. ljB5a,ll is identical to Eq. H2fia|l . so its solution is 
again inside the hard cores of the colloids, and given 
by Eq. I|28J1 outside the hard cores, provided that we fix 

^o = l/(l-»?). 

The second equation l|B5hJ) is not identical to its coun- 
terpart Eq. Il26hjl . The solution is quite similar though: 
we need to fix the hard core parameter fivo to 



(3v = -Z 
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in order to make sure that the charge density p(r) is a 
multi-centered sum of DLVO profiles. The solution (in 
k-space) is then given by 



p k =(27T) 3 {p+(p 
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which, in real space, is indeed a multi-centered sum 
p(r) = J2i Pi(l r — with the individual profiles given 
by Eq. Note that Eq. © and {EJ only differ in 

the oc fc 2 <5(k) term. As a consequence, the profiles p(r) 
resulting from those two equations, are identical, but the 
minimum of the functional differs. 

Upon insertion of the equilibrium density profiles into 
the functional {B2j, we immediately notice that the last 
term on the right hand side vanishes. The grand potential 
then becomes 
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The integral in this expression can now be calculated 
using Parseval's theorem and the expression <l A3|) for the 
Fourier transform of V(r). The result is 



(B4) _L|drp(r)/3y(r) = 




(B9) 
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so that the grand potential eventually becomes The volume term j3<& is exactly equal to the one that 

was found previously in Eqs. ffi7jl and iJHEJl; the colloidal 
( Ze Ra \ s--\ e ~ KR a p a i r interaction, however, reduces to a purely DLVO in- 

'"" "~ > 1 'Raj B teraction, i.e. the factor A we found before, is now equal 

(BIO) toO. 
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